function [lon,lat,ann_av,ann_std] = get_SETASAnnualdiag(modelvar)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% THIS FUNCTION COMPUTES THE  ANNUAL MEAN AND STANDARD DEVIATION OF THE 
% SELECTED SETAS MODEL VARIABLE.
%
% INPUT     : - modelvar  : a SETAS variable. Available options are 'temp or 
%                          'salt' 
% OUTPUT    : - lon       : SETAS longitude 
%             - lat       : SETAS latitude
%             - ann_av    : parameter annual mean
%             - ann_std   : parameter annual standard deviation
%
% Author: Benedicte Pasquer, IMOS/eMII (http://imos.org.au/)
% email: benedicte.pasquer@utas.edu.au
% May 2013   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

dirname = pwd ;
flist = dir(fullfile (dirname,'setas_simple_*.nc'));

% LOOP THROUGH SETAS FILE LIST TO EXTRACT SURFTACE VALUES)

for n = 1:length(flist)
% READ SETAS GRID AND TIME INFORMATION  FROM SETAS NETCDF MONTHLY FILE
    ncid = netcdf.open(fullfile(dirname,flist(n).name),'NC_NOWRITE');

    varid = netcdf.inqVarID(ncid,'time');
    time = netcdf.getVar(ncid,varid); 
    varid = netcdf.inqVarID(ncid,'latitude');
    lat = netcdf.getVar(ncid,varid); 
    varid = netcdf.inqVarID(ncid,'longitude');
    lon = netcdf.getVar(ncid,varid); 
    lon = lon';lat = lat';
% GET SPECIFIED VARIABLE DATA
    switch modelvar
        case 'temp'
            varid = netcdf.inqVarID(ncid,'temp');
            temp = netcdf.getVar(ncid,varid,'double');           
            var = temp(:,:,21:22,:);
            var = mean(var,3);
            var = squeeze(var);%REMOVE THE SINGLETON DIMENSION
            var = permute(var,[2 1 3]);
        case 'salt'
            varid = netcdf.inqVarID(ncid,'salt');
            salt = netcdf.getVar(ncid,varid,'double'); 
            var = salt(:,:,21:22,:);
            var = mean(var,3);
            var = squeeze(var); %REMOVE THE SINGLETON DIMENSION
            var = permute(var,[2 1 3]);

    end
    
    netcdf.close(ncid)
% STORE MONTHLY SURFACE DATA IN MATRIX;
    var_mthly_av(:,:,n) = mean(var,3);
end


var_mthly_av(var_mthly_av>1000) = NaN; 
%COMPUTE ANNUAL MEAN BASED ON THE 2YR DATASET. 
ann_av = mean(var_mthly_av,3);

%COMPUTE STD . NOTE THAT STD IS COMPUTED USING MONTHLY MEAN RESULTS NOT
%HOURLY RESULTS TO DISPLAY SEASONAL INSTAEAD OF HIGH FREQUENCY (hourly to
% daily) INFORMATION 
ann_std = std(var_mthly_av,0,3);





